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ABSTRACT 

The atmospheres of weakly-magnetized neutron stars expand hydrostaticaUy and 
rotate differentially during thermonuclear X-ray bursts. Differential rotation is proba- 
bly related to the frequency drifts of millisecond burst oscillations exhibited by about 
a dozen nuclear-powered X-ray pulsars. Here, we analyze the linear stability of this dif- 
ferential rotation with respect to local, axisymmetric, multi-diffusive MHD perturba- 
tions, at various heights in the neutron star atmosphere. Unstable magneto-rotational 
modes are identified from within to well above the burning layers. Properties of the 
fastest growing modes depend sensitively on the local magnetic field geometry. Linear 
estimates suggest that momentum transport due to magneto-rotational instabilities 
can affect atmospheric rotation profiles on time-scales relevant to burst oscillations. 
This transport would likely strengthen the coherence of burst oscillations and con- 
tribute to their drifts. 
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1 INTRODUCTION 

Since their discovery in the 1970s (see Lewin, van Paradijs & 
Taam 1995 for a review), thermonuclear X-ray bursts have 
been observed from a large number of low-mass X-ray bina- 
ries containing weakly-magnetized accreting neutron stars. 
The basic physical mechanism responsible for X-ray bursts, 
namely mass accumulation on the neutron star surface at 
rates such that a thermonuclear He or mixed H/He flash is 
eventually triggered, is well understood (see, e.g., Bildsten 
1998 and references therein). The phenomenology of X-ray 
bursts is rich, however, and theoretical models face some 
difficulties when detailed comparisons with the observations 
are attempted (e.g. Bildsten 2000; Galloway et al. 2003). 

This situation has become more serious, and more ex- 
citing, in the past decade, with the discovery of transient 
millisecond (300-600 Hz) oscillations during some X-ray 
bursts, the so-called burst oscillations (Strohmayer et al. 
1996). While the presence of oscillations during burst rise 
may be understood eis resulting from the rotational modula- 
tion of a nuclear burning spot spreading around the stellar 
surface (Spitkovsky, Levin & Ushomirsky 2002), it is un- 
clear why oscillations should be visible past the burst peak, 
when the burning fuel has presumably been entirely ignited. 
Also intriguing is the nature of frequency drifts exhibited 
by most burst oscillations (generally spin-ups, by a few Hz 



at most) and the origin of the diversity they show in terms 
of drift time-scales and amplitudes (see Strohmayer & Bild- 
sten 2003 for a review) . Observations of the burst oscillation 
phenomenon have been very rewarding. In particular, it has 
recently been shown that the burst oscillation asymptotic 
frequency (after complete spin-up) and the independently- 
known neutron star spin frequency are closely related in 
two millisecond accreting pulsars (Chakrabarty et al. 2003; 
Strohmayer et al. 2003). Still, our understanding of the phys- 
ical mechanisms involved in burst oscillations and associated 
frequency drifts remains limited. 

Following Strohmayer et al. (1997), Gumming & Bild- 
sten (2000) proposed that frequency drifts originate from 
the rotational evolution of heated burning layers, as they 
cool, hydrostaticaUy contract and nearly recover their orig- 
inal rotation, assuming their speciflc angular momentum is 
conserved at all times. These same authors later realized, 
however, that the largest observed frequency drifts are in 
excess of the amplitudes expected in their burning layer sce- 
nario, by as much as a factor 2-3. Apparently, only layers 
higher up in the atmosphere would have the right amount 
of rotational offset to explain the largest drift amplitudes 
(Gumming et al. 2002). More recently, Ghakrabarty et al. 
(2003) have also argued that the very rapid drift (spin-up) 
that they have observed during the rise phase of one of the 
bursts of SAX J1808-3658 effectively rules out the scenario of 
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Gumming & Bildsten (2000). There is currently no solid al- 
ternative explanation for the origin of burst oscillations and 
associated frequency drifts, even though both Rossby wave 
(Heyl 2003; but see Lee 2003) and zonal flow (Spitkovsky et 
al. 2002) scenarios have been proposed. 

In their study of neutron star atmospheres, Gumming 
& Bildsten (2000) have analyzed a number of mechanisms 
which could make the vertical difi^erential rotation estab- 
lished via burst-induced hydrostatic expansion unstable. 
They have not identified any clear destabilizing process 
(which would presumably lead to angular momentum trans- 
port between atmospheric shells) and this was the main mo- 
tivation for their assumption of conserved specific angular 
momentum. The purpose of the present study is to recon- 
sider the issue of stability of differential rotation in the at- 
mospheres of weakly-magnetized, rapidly-rotating neutron 
stars. 

As noted by Gumming & Bildsten (2000), the strong 
vertical thermal stratification of a neutron star atmosphere 
stabilizes it with respect to adiabatic perturbations, even in 
the presence of strong vertical differential rotation. Although 
not considered by these authors in their analysis, one in- 
deed verifies that the vertical differential rotation expected 
in the present context is stable according to both the hydro- 
dynamicaJ Solberg-H0iland adiabatic criteria (e.g. Tassoul 
1978) and their generalization for weakly-magnetized fluids 
(Balbus 1995). 

This situation is reminiscent of the solar radiative in- 
terior, where thermal stratification also plays a strongly 
stabilizing role. Goldreich & Schubert (1967) and Pricke 
(1968) have shown, however, that rotational instabilities can 
still exist, provided perturbed fiuid elements are allowed 
to exchange heat with their environment much faster than 
they exchange momentum. The stabilizing role of the atmo- 
spheric thermal stratification is then effectively neutralized 
when a perturbed fiuid element reaches thermal equilibrium 
with its environment while its original momentum remains 
largely unchanged. In the solar interior, radiative heat dif- 
fusion is indeed orders of magnitude faster than viscous dif- 
fusion of momentum, and this was the initial motivation 
for the double-diffusive analysis of Goldreich & Schubert 
and Fricke (hereafter GSF altogether). The basic double- 
diffusive mechanism invoked in these instabilities is not too 
different from that of "salt-finger" instabilities operating in 
the Earth's oceans, except that the destabilizing salinity 
stratification is replaced by a destabilizing angular momen- 
tum stratification in the stellar context. 

As we shall see below, conditions in the atmospheres 
of weakly-magnetized neutron stars during bursts are also 
strongly multi-diffusive. It is then possible that rotational in- 
stabilities will develop and affect the differentially-rotating 
atmospheric layers in such a way as to influence the ob- 
servational properties of burst oscillations. To address this 
question properly, it is necessary to account for the pres- 
ence of magnetic fields in neutron star atmospheres. Re- 
cently, Menou, Balbus & Spruit (2004) have generalized the 
work of GSF to weakly-magnetized fluids. They have shown 
that the combination of weak magnetic fields and a multi- 
diffusive situation may have contributed to establishing the 
state of near solid-body rotation inferred from helioseismol- 
ogy for the solar interior. 

In this paper, we present a speciflc application of 



the stability analysis of Menou et al. (2004) to weakly- 
magnetized, differentially-rotating neutron star atmospheres 
and we discuss possible consequences for the phenomenol- 
ogy of burst oscillations and associated frequency drifts. We 
describe our method of solution in §2. In §3, we present the 
results of our stability analysis at various heights in the neu- 
tron star atmosphere. The relevance of these results to burst 
oscillations is discussed in §4. 



2 METHOD 

Following Menou et al. (2004), we perform a local, lin- 
ear stability analysis of differential rotation in a weakly- 
magnetized, stably-stratified fiuid. In this section, we first 
recall the dispersion relation obeyed by axisymmetric modes 
in a general context in §2.1. We describe how this dispersion 
relation can be applied to the specific problem of weakly- 
magnetized neutron star atmospheres during bursts in §2.2. 
We then describe how we obtained numerical solutions to 
the dispersion relation in §2.3. In §2.4, we discuss the is- 
sue of the efficiency of momentum transport, that we esti- 
mate from linear mode properties. Our basic analysis follows 
closely that presented in Menou et al. (2004), where addi- 
tional details on the method can be found. 



2.1 Dispersion Relation 

We initially work in cylindrical coordinates, (R,(j),Z), but 
we will also later express results in terms of spherical coordi- 
nates, (r,6,(f)). We consider axisymmetric Eulerian pertur- 
bations with a WKB space-time dependence exp[i(fe ■ r — 
ijJt)], where k = (fcfl,0, fcz) = (fcr,fc9,0). Starting from the 
set of MHD equations (continuity, momentum, induction, 
energy) that include the effects of viscosity, resistivity and 
heat conduction 



dt 



+ V.(pi;) = 0, 



(1) 



— = V X (v X B)- 7?Vx (V X B) , 



1 „dlnPp- 



(7-1) dt 



(2) 
(3) 
(4) 



one derives, under the Boussinesq approximation and 
ignoring self-gravity, the following flfth-order dispersion re- 
lation for local, axisymmetric, multi-diffusive modes. 
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where 
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A developed form of the dispersion relation can be found 
in Menou et al. (2004). The Brunt- Vaisala frequency, A'^, 
which measures the magnitude of atmospheric thermal strat- 
ification, comes out of the first bracket term above. Similarly, 
the epicyclic frequency, which measures the angular momen- 
tum "stratification," comes out of the second bracket term. 
The notation is standard and identical to that adopted in 
Menou et al. (2004): v is the fluid velocity, p is the mass 
density, P is the pressure, "1> is the gravitational potential, 
T is the temperature, B is the magnetic field, fj. is the dy- 
namic viscosity, v — jj,/ p is the kinematic viscosity, rj is the 
resistivity, x is the heat conductivity and va is the Alfven 
speed. 

A value 7 = 5/3 is adopted for the gas adiabatic in- 
dex and a perfect gas equation of state is assumed. Mod- 
els of neutron star atmospheres during thermonuclear X-ray 
bursts show that radiation pressure sometimes contribute a 
non-negligible fraction to the total pressure, especially close 
to the burning layers (see, e.g.. Gumming & Bildsten 2000). 
We neglect this contribution here, but it should be included 
for better accuracy. Degeneracy is neglected, because it is 
lifted during bursts, and general relativistic effects are ig- 
nored. When discussing burning layers, we also ignore the 
effects that any remaining nuclear burning may have on sta- 
bility, both as a heat source and in modifying the composi- 
tion of burning regions. The basic state rotation is given by 
fl — {0,0,Q,{R, Z)) along the Z-axis. The basic state mag- 
netic field, whose geometry is specified below, is assumed to 
be weak with respect to both rotation and thermal pressure. 



2.2 Neutron Star Atmosphere During Burst 

In order to solve the above dispersion relation (Eq. [S|), we 
must specify values for all the physical parameters entering 
the problem. Often, we rely on the recent, one-dimensional, 
hydrostatic models of Gumming & Bildsten (2000) for neu- 
tron star atmospheres during bursts. We consider condi- 
tions in three representative layers, right after burst trig- 
ger and the hydrostatic expansion that follows it. The first 
such layer, labeled "Burning" in Table 1, is representative of 
conditions in the burning layers. The second layer, labeled 
"High", corresponds to typical conditions high in the atmo- 
sphere, well above the burning layers but still well below the 
stellar photosphere (conditions provided by A. Gumming, 
private communication). The third layer, labeled "Middle", 
corresponds to atmospheric conditions intermediate between 
the first two. Our choice of layers, although somewhat arbi- 
trary, will prove to be useful to illustrate variations of the 
stability properties with height in the neutron star atmo- 
sphere. We list in Table 1 the values of the density, p, tem- 
perature, T, and Rosseland-mean radiative opacity, k (from 
electron scattering), in the three layers of interest. 

We need to estimate the values of various microscopic 
diffusion coefficients relevant to the stability problem. From 



Spitzer (1962) and Schwarzschild (1958), assuming that 
the plasma is mostly composed of hydrogen, the (ion- 
dominated) dynamic viscosity, p, the radiative kinematic 
viscosity, i/r, the electric resistivity, i], the radiative heat 
conductivity, Xrad, and the corresponding heat diffusivity, 
Crad, are given by 



pu ~ 2.2 X 10"^^-;— r g cm~^ s~\ 
In A 



\rad 
^rad 



15 hip'^c^ ' 
5.2xl0"|^cm^ s-\ 

lerV 

c.g.s.. 



3kp 
7 - 1 T 



Xrad Cm^ S \ 



(6) 
(7) 
(8) 
(9) 
(10) 



where In A is the Goulomb logarithm, c is the speed of 
light and a is the Stefan-Boltzmann constant. We list in 
Table 1 the values of all three diffusivity coefficients for the 
three layers of interest. In addition, the values of the Prandtl 
number, = Vtot/£.rad = + i'r)/irad, and the "Acheson 
number," = rj/^rad, which together measure the hierar- 
chy of diffusive processes, are listed. Viscous diffusion clearly 
dominates over resistive diffusion in the neutron star atmo- 
sphere (a situation which is opposite to that in the Sun's 
radiative interior; Menou et al. 2004), but radiative heat 
diffusion is the strongest by many orders of magnitude (as 
is the case for the solar interior). It is this property of rapid 
heat diffusion that allows multi-diffusive modes to overcome 
the stabilizing inffuence of the atmospheric thermal strati- 
fication. Note also how the situation becomes more double- 
diffusive as one goes up in the atmosphere (since the values 
of and e,, decrease), and how radiative viscosity even- 
tually dominates over ion viscosity high enough above the 
burning layers. 

A number of global parameters for the magnetized, 
differentially-rotating neutron star atmosphere must be 
specified. We adopt a value Q, = 2x 10^ s~^ for the neutron 
star angular velocity, corresponding to a typical spin period 
of ~ 3 ms (spin frequency ~ 300 Hz). We adopt a value 
N/Q, ~ 10^ for the ratio of the Brunt- Vaisala frequency 
to the angular velocity (Gumming & Bildsten 2000). Note 
that we only include the thermal stratification in our def- 
inition of the Brunt- Vaisala frequency since the dispersion 
relation that we are solving was derived ignoring the effects 
of composition gradients. The issue of stabilization due to 
composition gradients is important and is further discussed 
in §4. 

We need to specify the level of differential rotation 
present in the atmosphere. After burst trigger, all the at- 
mospheric layers, from the burning regions to those located 
higher up, are heated up. This occurs on a thermal timescale, 
which is much in excess of a local sound crossing time (Gum- 
ming & Bildsten 2000) and, as a result, the layers expand 
hydrostatically. Assuming conservation of specific angular 
momentum during this phase, differential rotation between 
spherical shells should then be established at a level cor- 
responding to dXriQ./ dlnr = —2, in all three layers of in- 
terest. We further assume that no differential rotation is 
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Layer p T 

(g cm-3) (10** K) 



(cm^ g -*-) (cm^ 



-*-) (cm^ s ^) (cm^ s ^) 



(cm^ s-1) (Utot/^rad) iv/^rad) 



Burning lO'' 10 0.15 

Middle 3 X 10^ 5 0.17 

High 100 2 0.23 



88 34 0.13 

680 2.5 X 10^ 0.28 

2.1 X 10^ 4.4 X 10'' 1.1 



4.5 X 10^ 2.7 X 10""' 2.9 x lO"'' 

5.6 X 10'' 5.7 X 10-5 5 x 10"^ 
2.4 X 10^ 1.9 X 10-5 4.6 x 10"^° 



Table 1. Typical Atmospheric Conditions During Burst 



present within spherical shehs. It should be noted that any 
such negative differential rotation (i.e. dln^l/dO < 0) would 
be unstable to adiabatic magneto-rotational modes (Balbus 
1995; Menou et al. 2004). 

Finally, we must specify a geometry and a strength for 
the basic state magnetic field present in the atmospheric 
layers. Although other configurations are possible (for in- 
stance with high multipoles; see also Ruderman 2003), we 
will assume here, for simplicity, a misaligned dipole geom- 
etry. More precisely, because our analysis is axisymmetric 
in nature, we must use an "axisymmetric dipole" geometry, 
which best resembles that of a true dipole. The spherical- 
radial and polar components of the field are thus 

Br = Bp cos(6» - 6»„,i,), Bg = ^ sm{e - Sn,;,), (11) 

where 8 is the polar angle (^ = or vr on the neutron 
star rotation axis), 6'mis is the equivalent of the angle by 
which the dipolar magnetic axis is misaligned from the ro- 
tation axis and Bp is the polar field strength at the neutron 
star surface (assumed to be the same in all layers of inter- 
est). With this geometry, field lines have an axisymmetric 
configuration such that they are nearly vertical for angles 
9 ~ Smis (magnetic "pole," which is an axisymmetric band, 
really) and nearly horizontal at the magnetic equator. This 
setup allows us to study the effects of the magnetic field 
geometry on the stability properties. We have investigated 
stability in the three atmospheric layers of interest for values 
of Bp ranging from lO'^ to 10® G, considered as typical for 
weakly-magnetized accreting neutron stars. 



2.3 Numerical solutions 

The complexity of the dispersion relation (Eq. 0) makes 
it difficult to analytically derive necessary and sufficient 
conditions for stability (Menou et al. 2004). In addition, 
since we are interested in studying the properties of unstable 
modes, we want to solve the dispersion relation numerically, 
with the appropriate set of parameters. To do this, we use 
the Laguerre algorithm described by Press et al. (1992) to 
solve equation Q as a fifth-order complex polynomial for 
a = —iu), with real coefficients. 

We have found useful to rewrite both the rotational and 
thermal stratification bracket terms appearing in the disper- 
sion relation in spherical coordinates {r,9,(j)). The thermal 
stratification term then becomes function only of A'^^ , 9 and 
the radial and angular wavevectors, kr and kg, for a spheri- 
cally symmetric star. The rotational stratification term, on 



the other hand, explicitly depends on the amount of differ- 
ential rotation within and between spherical shells, that we 
express as din Q/ 89 and 91nf2/91nr, respectively. As we 
have mentioned earlier, we are interested here in the specific 
case corresponding to d\nQ/dd = and c?lnn/91nr = —2. 

For each of the three atmospheric layers of interest, 
we search for unstable modes (which correspond to a-roots 
with positive real parts) at various locations on the stellar 
sphere. We consider values of the polar angle, 9, in the range 
^ 7r/2 (pole to equator). At each location, we perform a 
search for unstable modes by varying the wavevectors kr and 
kg independently in the range ±[27r/Zmax, 27r/Zniin]. The first, 
large-scale limit, Zmax, is chosen to be ^ H, the pressure scale 
height, so that the local and weak field assumptions of our 
analysis remain valid. The second, small-scale limit, Zmin, is 
chosen to be 3> A, the mean free path, so that the analysis 
remains valid given our use of fluid equations. We estimate 
H ~ ?R.T/g, where 5R is the perfect gas constant, A ~ u/VWT 
and we use (/ ~ 2 x 10^* cm s"^ as a typical neutron star sur- 
face gravity. In addition, we systematically perform focused 
searches in regions of the wavevector space corresponding to 
nearly cylindrical-radial displacements (|fefl/fcz| — > 0) and 
nearly vertical displacements Qku/kzl oo). 



2.4 Efficiency of Transport 

The differential rotation established in neutron star at- 
mospheres during bursts is an intrinsically transient phe- 
nomenon. Absent any transport mechanism, the various 
differentially-rotating spherical shells would essentially re- 
cover their original rotation (modulo a small offset due to 
changes in molecular weight; Gumming & Bildsten 2000) 
once the atmosphere cools down to pre-burst conditions, 
typically after several tens of seconds. In the presence of 
a transport mechanism, however, the issue of the efficiency 
of momentum transport is of direct relevance to the observa- 
tions because a mechanism effective enough to act on time- 
scales of tens of seconds or shorter would influence the ro- 
tational evolution of atmospheric layers during burst. 

It is unclear whether one can determine, from a linear 
analysis only, the efficiency of transport that would result 
from non-linear growth and turbulence. We will adopt an 
often-used prescription here, by assuming that the efficiency 
of transport can be crudely estimated from the growth rate 
and length-scale of the fastest-growing unstable linear mode. 
There is no general justification for this simple prescription, 
except perhaps examples of physical situations for which 
it is not too far off (for instance, the emergence of mid- 
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latitude weather structures in the Earth's turbulent atmo- 
sphere with sizes matching linear predictions for the fastest 
growing baroclinic modes is considered a success of baro- 
clinic theory; see Pedlosky 1987). Our goal here is not to 
accurately estimate the efficiency of momentum transport 
between shells, but merely to question the observational rel- 
evance of this transport given typical burst durations. 

It should be noted that the non-linear numerical simu- 
lations of Korycansky (1991) lend support to the simple pre- 
scription we adopted. This author found, for an equatorial 
setup (i.e. aligned thermal and angular momentum strati- 
fications), that the non- linear behavior of hydrodynamical 
double-diffusive (GSF) modes is well described by a com- 
bination of linear growth rate and spherical-radial length 
scale of the instability. Assuming this property carries over 
to magnetized modes and to other geometries (i.e. away from 
the stellar equator), we will be estimating the efficiency of 
momentum transport across spherical shells with the quan- 
tity crmax/(A:r)max (dimeusiou of a viscosity), where (Tmax is 
the growth rate of the fastest-growing mode and (fcr)max 
is the spherical-radial wavevector component of that same 
mode. Clearly, using linear theory to predict a non-linear 
outcome is a daring extrapolation, so one should treat our 
results on transport efficiencies with much caution. 



3 RESULTS 

For each of the three atmospheric layers of interest, we have 
run several models to identify parameter dependencies. To 
best illustrate these dependencies, we have chosen a specific 
subset of models, listed in Table 2. In most cases, only the 
magnetic field strength. Bp, or the value of Omis, the angle by 
which the magnetic "pole" is misaligned from the rotation 
axis, are varied. But we also show the effects of variations 
in the neutron star rotation rate, Jl, the rate of spherical- 
radial differential rotation, 91nn/91nr, and the maximum 
allowed scale for perturbations. 

Quite generally, many unstable modes were found 
in the searched wavevector space, as expected for the 
strongly multi-diffusive conditions considered. All the un- 
stable modes that we have identified were of the "direct" 
type (i.e. with Im{a) = 0), as opposed to the other multi- 
diffusive "overstable" type (when \Im{a)\ > Re{a) > 0). 
This is a consequence of the strongly stabilizing thermal 
stratification and of the fast rate of heat diffusion (relative to 
momentum diffusion), as was the case for the solar radiative 
zone (Menou et al. 2004). To simplify the description of these 
unstable modes, we focus our attention on the properties of 
the fastest-growing mode identified, as a function of location 
on the stellar sphere (i.e. polar angle, 0). While we cannot 
guarantee that we have identified the single fastest-growing 
mode at any location, convergence tests with increasingly 
higher resolutions for the wavevector grid and the polar an- 
gle grid were successfully passed, suggesting that our search 
method is robust. It is worth pointing out that we have often 
experienced unsatisfactory numerical convergence for small 
values of the polar angle, i.e. ^ ~ (close to the rotation 
axis). This may be caused by the condition for a weak mag- 
netic field compared to rotation becoming increasingly more 
stringent at small 8 values (where pressure gradients become 
nearly orthogonal to the centrifugal force), but we have not 
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Figure 1. Growth rate of the fastest-growing mode, ffmax, in 
units of angular velocity, Q, as a function of polar angle, Q (rota- 
tional equator at 6 = 7r/2). The solid, dotted, short-dashed and 
long-dashed lines correspond to models B3, Bl, B2 and HI in 
Table 2, respectively. 

explored this possibility further. We also do not explicitly 
show results for a magnetic field strength Bp — 10^ G. For 
reasons that will become clear by the end of this section, 
searches for unstable modes in this stronger field situation 
have been unsuccessful high in the atmosphere. 

We begin with a discussion of the general properties 
of fastest-growing unstable modes identified in the burning 
layers, before comparing them with the properties of fastest- 
growing modes identified higher-up in the atmosphere. 

3.1 Burning Layers 

Figure shows the maximum growth rate among all un- 
stable modes, in units of the neutron star angular velocity, 
f2, as a function of polar angle, Q. The solid, dotted and 
short-dashed lines correspond to models B3, Bl and B2 in 
Table 2, respectively. These models differ in the strength of 
the magnetic field and by how much the magnetic "pole" is 
misaligned from the rotation axis at = 0. 

In all three models, the growth rate of the fastest- 
growing mode is a function of polar angle, d, which shows 
that the combination of magnetic field and stellar spheri- 
cal geometries have an infiuence on stability properties. The 
similarity between the dotted and short-dashed lines also 
indicates that the magnetic field strength does not strongly 
affect growth rates, although we will see later that it does 
matter in a global sense (§3.2). On the other hand, there is 
a clear distinction between the solid line (model B3) and the 
dotted and short-dashed lines, indicating that the magnetic 
field geometry alone infiuences the stability properties quite 
a lot. 

These trends are confirmed by looking at the wavenum- 
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Model Layer Bp 6^;^ Special 

(Gauss) (degrees) 



Bl 


Burning 


-2 


108 







B2 


Burning 


-2 


10^ 







B3 


Burning 


-2 


108 


30 




Ml 


Middle 


-2 


108 


30 




M2 


Middle 


-2 


108 







M3 


Middle 


-2 


108 


20 




M4 


Middle 


-2 


108 


30 


n X 2 


HI 


High 


-2 


10^ 


30 




H2 


High 


-1 


10^ 


30 




H3 


High 


-2 


10^ 


30 


large imax 



Table 2. Models and Parameters 
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[radians] 



Figure 2. Ratio of spherical-radial to polar wavenumber, \kr/kg\, 
for the fastest-growing mode, as a function of polar angle, 6. No- 
tation is identical to Fig. 1, for the same four models (solid: B3; 
dotted: Bl; short-dashed: B2; long-dashed: HI). 

ber properties of the same modes. Figure |5| shows the ra- 
tio of spherical-radial to polar wavenumber, \kr/ke\, for the 
fastest-growing modes, as a function of polar angle, 9. The 
notation is identical to Figure for ease of comparison. 
Once again, the magnetic field strength is found to have lit- 
tle influence on the mode properties, but the magnetic field 
and stellar spherical geometries do.^ 

^ We have checked that the discontinuities present in the short- 
and long-dashed lines in figure|5] and absent from the correspond- 
ing growth rates in figure^ persist even at higher wavenumber or 



Figure |5| reveals an interesting wavenumber property 
of fastest-growing modes. In the two models with a mag- 
netic dipole aligned with the rotation axis, the nature of the 
modes changes from mostly spherical-radial {\kr/ke\ <C 1) 
for small 6 values to mostly horizontal (|fcr/fce| 3> 1) as 
one approaches the stellar equator {6 it/ 2). In the model 
with a magnetic "pole" misaligned by Smis ~ 30 deg, on the 
other hand , \kr/ke\ is ^ 1 on a fair fraction of the stellar 
sphere, although it approaches unity at the stellar equator. 
This behavior can be understood as resulting mostly from 
the magnetic field geometry. At the magnetic "poles," where 
the field is nearly spherical-radial, the fastest-growing mode 
tends to be nearly spherical-radial itself,^ with |fcr/fee| ^ 1. 
As one approaches the magnetic equator, however, the field 
becomes nearly horizontal and so does the fastest-growing 
mode. This is consistent with a tendency for the fastest- 
growing mode to minimize the amount of magnetic tension 
it is subject to, by favoring displacements nearly along field 
lines. 

Figure|3]shows estimates of the efficiency of momentum 
transport in the spherical-radial direction, based on fastest- 
growing linear mode properties, as a function of polar angle 
6. Again, notation is similar to that in Figs. 1 and 2 for 
models Bl (dotted line), B2 (dashed line) and B3 (solid 
line). In model B3, which corresponds to a misaligned mag- 
netic dipole, the estimated efficiency of transport is strongly 
peaked at the magnetic "pole" (because of a significant in- 
crease of the mode radial wavelength, as we shall see below 
and in Fig. The horizontal dash-dotted line in Figure |H] 
indicates the efficiency typically required in the burning lay- 
ers to affect the rotation profile over a local pressure scale- 
height, H ~ ^T/g, in 1 second. This reference efficiency 
was simply calculated as cm^ s~^. It would not make 

polar angle grid resolutions. Consequently, they appear to have 
a physical origin (possibly a change in the nature of the fastest- 
growing mode with local geometry) . We have not investigated this 
further. 

^ The nearly incompressible modes described by our dispersion 
relation are characterized by displacements orthogonal to the 
mode wavevector, fc. 
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Figure 3. Estimate of the efficiency of spherical-radial momen- 
tum transport, (Tmax/k^ (dimension of a viscosity), based on 
fastest-growing mode properties in typical burning layer condi- 
tions, as a function of polar angle, 9. The solid, dotted and dashed 
lines correspond to models B3, Bl and B2 in Table 2, respectively. 
The location of the magnetic "pole" at = 30 deg in model B3 is 
obvious. The horizontal dash-dotted line indicates the efficiency 
required to affect the rotation profile over a local scale-height in 
1 s. 



much sense to calculate a reference efficiency for a length- 
scale exceeding H because both atmospheric and stability 
conditions start to change significantly over a local pressure 
scale-height. 

Interestingly, the maximum momentum transport effi- 
ciency estimated from fastest-growing mode properties, at 
the magnetic "pole," falls short of the reference efficiency 
only by a factor of a few, suggesting that several seconds 
may be sufficient for the rotational profiles to be affected. 
Away from the magnetic "pole" (and essentially everywhere 
on the stellar sphere in aligned-dipole models Bl and B2), 
estimated transport efficiencies appear too small, by several 
orders of magnitude, to affect rotation profiles over a local 
scale-height on observationally relevant time-scales. In fact, 
the efficiencies indicated by the dashed and dotted lines in 
Figure 121 are small enough to be comparable to that of mi- 
croscopic viscosity {i/ + u,.) in the burning layers, as can be 
seen from Table 2. 

Many of the basic properties we have just described for 
fastest-growing modes in burning layers carry over to higher 
atmospheric layers. There are important differences, how- 
ever, that we will now illustrate. Rather than going through 
a description of growth rates and wavenumbers indepen- 
dently for all cases of interest, we will often summarize re- 
sults by focusing on transport efficiencies (which are a com- 
bination of the two) and how these are affected by changes 
in physical conditions higher up in the atmosphere. 



10 I 1 < < r 
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Figure 4. Estimate of the efficiency of spherical-radial momen- 
tum transport for typical conditions in a high-atmospheric layer, 
as a function of polar angle, 9. The solid, dotted and dashed lines 
correspond to models Hf, H3 and H2 in Table 2, respectively. 
The horizontal dash-dotted line indicates the efficiency required 
to affect the rotation profile over a local scale-height in f s. 



3.2 High Atmosphere 

In addition to the three burning layer models we have al- 
ready discussed, figures and |5| show, as long-dashed lines, 
the growth rates and wavenumber properties of fastest- 
growing modes in our fiducial model for a high atmospheric 
layer (model HI in Table 2). Although strong similarities 
exist between models HI and Bl (solid line), especially in 
terms of wavenumber properties, a striking difference is the 
absence of any unstable mode in the vicinity of the mag- 
netic "pole," at ^ ~ 30 deg. As we shall see, this results 
from the length-scale of potentially unstable modes exceed- 
ing the maximum value allowed by our search, Zmax ~ H 
(the local scale height). 

Figure |1] shows estimates of the efficiency of momen- 
tum transport in the spherical-radial direction, based on 
fastest-growing linear mode properties in the high atmo- 
spheric layer, as a function of polar angle 0. The solid, dotted 
and dashed lines correspond to models HI, H3 and H2, re- 
spectively. The horizontal dash-dotted line indicates again 
the efficiency typically required to affect the rotation profile 
over a local scale height in 1 s (~ cm^ s~^). 

Model H3 (dotted line) is voluntarily artificial, in that 
it allows the length-scale of unstable modes to be as large 
as the neutron star radius, _Rns (taken to be 10 km) , rather 
than the local scale-height of the high-atmospheric layer, 
H ~ 150 cm. The efficiency of momentum transport in 
this artificial model is strongly peaked around the magnetic 
"pole," at S ~ 30 deg, and much resembles that of model 
B3 (solid line in Fig. |3 . Clearly then, the difference with 
model HI, and in particular the absence of any unstable 
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mode in the vicinity of the magnetic "pole," results from 
the tendency of these modes to have large length-scales, in 
excess of /max ~ H, for typical conditions high-up in the 
atmosphere. 

This property is consistent with the gradual breakdown, 
as one goes up in the atmosphere, of the important weak- 
field assumption made in our analysis. Indeed, while the 
magnetic field strength may not vary significantly over a 
few scale heights, density and pressure do. As a result, our 
working hypothesis of a magnetic field weak compared to 
thermal pressure becomes increasingly difficult to satisfy. 
Specifically, while the Alfven speed, va oc Bp~^^^, increases 
in the more tenuous high atmosphere, the sound speed, 
cs oc T^^^, decreases and the requirement va ^ cs becomes 
less and less valid. It is only marginally satisfied in the high- 
atmospheric layer we have considered (va/cs ~ 0.1 for the 
values of Bp adopted) and this is the origin of the disap- 
pearance of unstable modes around the magnetic "pole" in 
model HI (e.g. Fig.^J- While large-scale displacements sub- 
ject to weak enough magnetic tension for instability could 
be found in principle, they are clearly outside the regime 
of validity of our weak-field analysis. Even higher up in the 
atmosphere, the weak field regime breaks down completely 
and all unstable diffusive magneto-rotational modes disap- 
pear. 

Still, some of the fastest-growing modes identified in 
our reahstic high-atmosphere HI model (solid line in Fig.^J 
could affect rotation profiles over a local scale height on 
rather short time-scales, especially near the stellar equa- 
tor. Provided that these unstable modes act to reduce the 
amount of differential rotation initially present, a rotational 
recoupling of unstable layers in the high atmosphere would 
then be expected on observationally interesting time-scales. 
It is important to notice, however, that the magnitude of mi- 
croscopic viscosity {v + Vr) in high-atmospheric layers (com- 
pare Table 1 and Fig. 2J is sufficient by itself to reduce the 
level of differential rotation on similar time-scales, so that 
the existence of unstable modes, although important, is not 
required to achieve rotational recoupling in these high lay- 
ers. 

Model H2 (dashed line in Fig. |1J , differing from model 
HI only by a smaller rate of spherical-radial differential ro- 
tation, d\nQ./d\iir — —1, shows that our conclusions on 
the stability and transport properties of diffusive magneto- 
rotational modes do not strongly depend on the exact level 
of differential rotation assumed to be present. 

Finally, for completeness, let us mention that, by en- 
forcing Bp = 0, we have been able to identify purely hydro- 
dynamical multi-diffusive (GSF) modes in high-atmospheric 
layers. These modes were characterized by systematically 
lower momentum transport efficiencies than their magne- 
tized counterparts. They did not exist in lower atmospheric 
layers, probably because local conditions are less double- 
diffusive as one approaches the burning layers (see Table 1). 
This result is consistent with the findings of Menou et al. 
(2004): a weak magnetic field acts as a catalyst for rotational 
instability in a multi-diffusive context, much like it does in 
an adiabatic context. 
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Figure 5. Estimate of the efficiency of spherical-radial momen- 
tum transport for typical conditions in a middle-atmospheric 
layer, as a function of polar angle, d. The solid, dotted, long- 
dashed and short-dashed lines correspond to models Ml, M4, M2 
and M3 in Table 2, respectively. The horizontal dash-dotted line 
indicates the efficiency required to affect the rotation profile over 
a local scale-height in 1 s. 
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Figure 6. Ratio of spherical-radial wavelength to local scale 
height, noted \kii/kr\, for the fastest-growing modes in a mid- 
dle atmospheric layer, as a function of polar angle, 0. Notation 
is identical to Fig. 5, for the same four models (Ml: solid; M2: 
long-dashed; M3: short-dashed; M4: dotted). 
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3.3 Middle Atmosphere 

Having identified important differences between unstable 
modes in burning and high-atmospheric layers, it is only 
natural to wonder what would the properties of these modes 
be in intermediate middle-atmosphere conditions. Figure |K| 
shows estimated efficiencies of momentum transport, based 
on fastest-growing linear mode properties, in four different 
models for middle-atmospheric conditions. The solid, dot- 
ted, long-dashed and short-dashed lines correspond to mod- 
els Ml, M4, M3 and M2 in Table 2, respectively. The hori- 
zontal dash-dotted line shows the typical efficiency required 
to affect rotation profiles over a local scale height in 1 s. Fig- 
ure El shows, for the same four models (and same notation 
as in figure |^, the ratio of the spherical-radial wavelength 
of the fastest-growing mode to the local scale height, noted 

\kH/kr\. 

Results for middle-atmospheric layers shown in figures|3 
and 1^ are indeed consistent with the conclusions reached 
for burning and high-atmospheric layers. In all three mod- 
els with a misaligned dipole (Ml, M3 and M4; Table 2), 
the effect of the magnetic field geometry is readily seen. 
In the vicinity of the magnetic "pole," fastest-growing dif- 
fusive magneto-rotational modes are characterized by pref- 
erentially spherical-radial displacements (nearly along field 
lines), with a radial wavelength reaching the maximum 
value, /max ~ H (Fig. EJ, because the weak magnetic field 
assumption is not ideally satisfied even for these middle- 
atmospheric conditions. As a result, the modes left unstable 
in the vicinity of the magnetic pole, which marginally sat- 
isfy the constraint on /max, have slower growth rates and re- 
duced efficiencies of momentum transport. This is the origin 
of the dips in figure |^ and plateaus in figure El at ^ ~ Omis- 
The model with an aligned dipole, on the other hand, is 
characterized by fastest-growing modes with smaller radial 
wavelengths and transport efficiencies. 

What makes middle-atmospheric layers interesting is 
that, in the models with a misaligned magnetic dipole, es- 
timated efficiencies of momentum transport are in excess of 
the reference value indicated by the horizontal dash-dotted 
line, in the vicinity of ~ 6'mis. These estimates are also 
much in excess of the transport caused by microscopic vis- 
cosity in these layers [utot = v + Ur 3 x cm^ s~^; 
Table 1). Consequently, the influence of these diffusive 
magneto-rotational modes is expected to dominate the ro- 
tational evolution of middle-atmospheric layers, on short 
enough timescales to possibly have observational conse- 
quences. This conclusion is not very sensitive to the exact 
value of the neutron star angular velocity, Q. A compari- 
son between models Ml and M4 (solid and dotted lines in 
Fig. EJ shows that a faster stellar rotation leads to faster 
growth and transport, but only moderately so. 



4 DISCUSSION AND CONCLUSION 

The atmospheres of weakly-magnetized neutron stars are 
subject to vertical differential rotation during thermonu- 
clear X-ray bursts. Conditions in these atmospheres are also 
strongly multi-diffusive and thus favor the development of 
magneto-rotational instabilities. We have presented a lin- 
ear stability analysis of the expected differential rotation. 



with respect to local axisymmetric MHD perturbations, in 
three representative atmospheric layers. Unstable diffusive 
magneto-rotational modes were identified in much of the 
differentially-rotating atmosphere. This indicates that mo- 
mentum transport between the various atmospheric layers 
could occur. 

Our main motivation to perform a stability analysis is 
the possibility that the transport resulting from instabilities 
may have observational consequences for millisecond burst 
oscillations and associated frequency drifts seen so far in a 
dozen accreting neutron star systems (nuclear-powered X- 
ray pulsars; e.g. Strohmayer & Bildsten 2003). The differen- 
tial rotation established during thermonuclear X-ray bursts 
is transient, lasting perhaps several tens of seconds (a typi- 
cal burst duration) . Therefore, it is not sufficient to identify 
magneto-rotational instabilities in this context. For these in- 
stabilities to have any effect, they must act to modify atmo- 
spheric rotation profiles on short enough timescales. We have 
attempted to address this issue by using linear estimates of 
the efficiency of momentum transport and have concluded 
that, indeed, these instabilities could be efficient enough to 
be of relevance to the bust oscillation phenomenon. 

There are reasons to believe that diffusive magneto- 
rotational instabilities will, in a stellar context, tend to re- 
duce the level of existing differential rotation and perhaps 
make the system approach a state of solid body rotation 
(see discussion in Menou et al. 2004; see also Korycansky 
1991). If correct, this conjecture would have interesting con- 
sequences for the coherence of burst oscillations. Indeed, it 
was noted by Cumming & Bildsten (2000) that the strong 
coherence of some burst oscillations implies, in the absence 
of any momentum transport between atmospheric shells, 
that the feature at the origin of the oscillations must be 
confined to a surprisingly thin layer (much thinner than a 
local scale height in their burning layer scenario). This fea- 
ture would otherwise be smeared out by differential rotation 
and the oscillations would lose their coherence. 

If, however, magneto-rotational instabilities act to re- 
duce the amount of differential rotation present during 
bursts, as we have suggested, they would effectively pro- 
mote the coherence of burst oscillations. Constraints on the 
thickness of the layer from which the modulated signal orig- 
inates would then be relaxed. Magneto-rotational instabili- 
ties would also contribute to the observed frequency drifts of 
burst oscillations if they act efficiently enough to modify the 
rotational evolution of atmospheric layers. This possibility 
is consistent with the results of Cumming et al. (2002) , who 
concluded that conservation of specific angular momentum 
cannot account for the largest observed frequency drifts. 

As attractive as it is, this scenario suffers from a num- 
ber of significant uncertainties. For instance, we have cau- 
tioned about the risks of estimating efficiencies of momen- 
tum transport from linear mode properties only. Also, the 
feature which is at the origin of burst oscillations must be 
non-axisymmetric. The assumption of axisymmetry made 
in our analysis could therefore be a limitation. We have 
illustrated how the properties of fastest growing magneto- 
rotational modes strongly depend on the local magnetic field 
geometry. This may be viewed as a model uncertainty be- 
cause the global field topology on a given neutron star is 
a priori unknown. On the other hand, this sensitivity to 
magnetic field geometry might in principle be related to the 
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link recently established between the phases of persistent 
and burst oscillations, because it suggests a special role for 
magnetic polar caps in the burst oscillation phenomenon 
(Chakrabarty et al. 2003; Strohmayer et al. 2003). One then 
wonders whether the diversity in drift timescales could be 
attributed to different field geometries in different systems, 
perhaps caused by varying degrees of field burying under 
continuous accretion (Gumming, Zweibel & Bildsten 2001). 

Another possible limitation of our work comes from ne- 
glecting the stabilizing role of composition gradients. In the 
vicinity of burning layers, substantial composition gradients 
may be expected in ashes, following the main episode of nu- 
clear burning. The stratification due to composition gradi- 
ents is important because, contrary to thermal stratification, 
it is not sensitive to a fast rate of heat transfer. Although 
it is difficult to assess the importance of stabilization due 
to composition gradients without a detailed study, the spe- 
cific example described by Goldreich & Schubert (1967), for 
a purely hydrodynamical situation, suggests that this effect 
should be included when discussing neutron star atmosphere 
during bursts. Another motivation to elucidate the role of 
stabilizing composition gradients is the possibility that non- 
axisymmetric versions of diffusive magneto-rotational modes 
could contribute to the friction term entering the burning 
front models of Spitkovsky et al. (2002). 

Rather than focusing on burning layers only, we have 
performed stability analyses at various heights in the 
differentially-rotating neutron star atmosphere. One moti- 
vation for this is the stabilizing effect of composition gradi- 
ents that we have just discussed. Even if such stabilization 
occurs in the vicinity of burning layers, magneto-rotational 
instabilities may still be able to operate higher-up in the at- 
mosphere. A second motivation is the indication, from the 
work of Gumming et al. (2002; see, e.g., their figure 4), that 
layers above the burning regions are involved in the burst 
oscillation phenomenon, since they are the only ones having 
the right amount of rotational offset to explain the largest 
observed frequency drifts. 

Our analysis for high atmospheric layers shows that 
magneto-rotational instabilities cease gradually to exist as 
one approaches the stellar photosphere, because the weak 
magnetic tension they require becomes increasingly difficult 
to satisfy. Still, according to our middle-atmospheric models, 
there is a region above the burning layers in which magneto- 
rotational instabilities can operate and possibly lead to a 
rather efficient transport of momentum. While burst oscil- 
lations have traditionally been associated with burning lay- 
ers, it may be worth giving more attention to these middle- 
atmospheric layers. Not only do they have the right amount 
of rotational offset to explain large frequency drifts, but 
these layers are also more advantageous for the visibility 
of oscillations, with a short heat diffusion time compared to 
the time they take to revolve around the star (see Gumming 
& Bildsten 2000 for a discussion of the visibility of burst 
oscillations) . 

We note that Gumming & Bildsten (2000) and Gum- 
ming et al. (2002) have invoked the mechanism of mag- 
netic field wind-up described by Spruit (1999) as a possible 
source of recoupling between differentially-rotating layers. 
In regions with weak enough stratification from composi- 
tion gradients for magneto-rotational instabilities to oper- 
ate, this process may not be relevant. Indeed, in the case 



of differentially-rotating accretion disks, it has been argued 
analytically and demonstrated numerically that the linear 
growth of the azimuthal field component caused by differ- 
ential rotation is not important because of the exponen- 
tial growth of magneto-rotational modes (Balbus & Ifawley 
1991; 1998). By analogy, one may expect the same linear 
growth to be irrelevant in the stellar context, given the exis- 
tence of exponentially-growing diffusive magneto-rotational 
modes. 

Beyond simple qualitative points, it is difficult to 
deepen an interpretation of burst oscillations and associated 
frequency drifts in terms of a rotational evolution which is 
partly driven by magneto-rotational instabilities. A time- 
dependent model for the coupled thermal, chemical and ro- 
tational evolution of a neutron star atmosphere during burst 
seems to be required to make more definite predictions. Such 
time-dependent models, ignoring the rotational evolution, 
already exist and exhibit complex behaviors. For example, 
transient convection early in the burst leads to some mixing 
of momentum and elements, while inverted composition gra- 
dients (which would presumably promote instabilities) are 
sometimes obtained after the main nuclear burning phase 
(see, e.g., Woosley & Weaver 1984; Woosley et al. 2003). 
In order to include the rotational evolution, these models 
would require reliable estimates of the efficiency of momen- 
tum transport (and of the resulting chemical mixing) , which 
could presumably be obtained from fully turbulent numer- 
ical simulations of the non-linear development of diffusive 
magneto-rotational instabilities. Until such advanced tools 
are developed, it is likely that observations will remain our 
best guide for understanding burst oscillations and associ- 
ated frequency drifts. 
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